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The recently developed density matrix quantum Monte Carlo (DMQMC) algorithm stochastically samples 
the A'^-body thermal density matrix and hence provides access to exact properties of many-particle quantum 
systems at arbitrary temperatures. We demonstrate that moving to the interaction picture provides substan¬ 
tial benefits when applying DMQMC to interacting fermions. In this first study, we focus on a system of 
much recent interest: the uniform electron gas in the warm dense regime. The basis set incompleteness error 
at finite temperature is investigated and extrapolated via a simple Monte Carlo sampling procedure. Finally, 
we provide benchmark calculations for a four-electron system, comparing our results to previous work where 
possible. 


I. INTRODUCTION 

The overwhelming majority of electronic structure 
studies of matter have been conducted at zero tempera¬ 
ture. This state of affairs has been justified as typically 
one is interested in the low-energy properties of con¬ 
densed matter systems or the room-temperature prop¬ 
erties of chemical systems. Due to recent experimental 
advances, however, there has been renewed interest in 
the thermodynamic properties of electron plasmas such 
as those found on the pathways to inertial confinement 
fusion^, in the interiors of Jupiter and other gas giants^, 
at the surfaces of solids after laser irradiation^, and in 
plasmonic catalysis"^. 

Of fundamental importance to the theoretical under¬ 
standing of these systems is the uniform electron gas 
(UEG), which has been pivotal to the development of 
modern quantum mechanical approaches to the low- 
temperature chemistry and physics of molecules and 
solids®’®. At finite temperatures the UEG can be de¬ 
scribed in terms of the density parameter, r^, and the 
degeneracy temperature, 0 = T/Tp, where Tp is the 
Fermi temperature. When both and 0 are of order 
one the system is said to be in the warm dense regime, 
with quantum, thermal and interaction effects all being 
important. It is in this region where analytical meth¬ 
ods, such as Ref. 7, are least useful and computational 
approaches such as quantum Monte Garlo are most im¬ 
portant. 

In a pioneering study. Brown et al.^ provided the first 
accurate data for the UEG in the warm dense regime us¬ 
ing the restricted path integral Monte Garlo (RPIMC) 
method®, from which the first accurate parameteriza- 
tions of finite-temperature density functionals have been 
produced^®’^®. Recent configuration path integral Monte 
Carlo (CPIMC) results^® have called into question the va¬ 
lidity of the restricted path approximation used in Ref. 


8, with significant disagreement between the two meth¬ 
ods at high densities and low temperatures. Simulations 
using a third method such as DMQMC would help to 
resolve this discrepancy^®. 

The exact thermodynamic properties of the UEG can 
be determined from the (unnormalized) iV-particle den¬ 
sity matrix 


where /3 = l/kpT. A direct evaluation of Eq. (1) requires 
all eigenvalues and eigenstates of H to be known and is 
an impossible task in all but the simplest of systems. 
The infinite basis set limit can be approached by only 
including determinants that can be constructed using a 
finite basis set of M plane waves, reducing the problem 
to the diagonalization of an x matrix. Even this 
problem is only tractable for very small M and N. 

Recently, Booth et al. have shown, through the de¬ 
velopment of the full configuration interaction QMC 
(FCIQMC) method, that full configuration interaction 
(FCI) quality results can be obtained at zero tempera¬ 
ture with no prior knowledge of the nodal structure of 
the wavefunction^"^. FCIQMC also often offers a substan¬ 
tially reduced memory cost compared to conventional 
FCI calculations. This has most dramatically been seen 
in the case of the UEG, where a space of 0(10^®®) Slater 
determinants was successfully sampled using the initia¬ 
tor adaptation of FCIQMC^®’®®. Subsequently, three of 
us used these ideas to develop the finite-temperature ana¬ 
logue of FCIQMC: density matrix quantum Monte Carlo 
(DMQMC). This was applied to the Heisenberg model as 
a proof of concept^^, but has not previously been used to 
study more realistic systems. 

Here we show how DMQMC can be applied to 
fermionic systems, starting with the UEG, thus opening 
the door to providing accurate, unbiased thermodynamic 
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results for problems of chemical interest. We note that 
CPIMC^® and Krylov-projected FCIQMC^® will likely be 
complementary approaches to both DMQMC and PIMC 
in the treatment of real systems. 

In Section II we outline the DMQMC method and 
show how moving to the interaction picture provides 
substantial improvements in statistical accuracy when 
treating weakly-correlated systems. In Section III we 
discuss basis-set extrapolation at finite temperatures in 
detail. In Section IV we present benchmark results 
for a four-electron system across the relevant parameter 
space, comparing, where possible, to previous results. Fi¬ 
nally, in Section V, we outline the limitations and future 
prospects of DMQMC. 

II. DENSITY MATRIX QUANTUM MONTE CARLO 

In this section we briefly outline the DMQMC algo¬ 
rithm; a more complete description is available in Ref. 17. 
DMQMC is applicable to any Hamiltonian but here we 
focus on the specific example of the UEG. We then ex¬ 
plain how to sample the density matrix in the interac¬ 
tion picture, show that this overcomes sampling issues 
found when treating weakly-correlated systems, and in¬ 
troduce a simple Monte Carlo scheme for sampling non¬ 
interacting density matrices in the canonical ensemble. 
Hartree atomic units are used throughout. 

A. Theory 

The unnormalized density matrix in Eq. (1) obeys the 
symmetrized Bloch equation 

% = ( 2 ) 
which can be solved using a simple Euler update scheme: 

p(/3 + A/3) = p(/3) - ^{um + mH) + 0(A/32). 

( 3 ) 

In DMQMC, Eq. (3) is solved stochastically by evolving 
a population of signed ‘psi-particles’, or ‘psips’, in a dis¬ 
crete operator space made of tensor products of Slater 
determinants. To this end, we rewrite Eq. (3) in matrix 
form: 

P\0 + A/3) = pij(/3) - SSik)pkj- (4) 

^ Pik(77kj - 5'^kj)] 

= Pij(/3) 3 ^ ^^(TfkPkj + Pik^kj)) (5) 

k 

where py = {Di\p\Dj), \Di) is a Slater determinant in 
the finite but large basis set. S' is a variable shift intro¬ 
duced to control the psip population^^’^°, and the last 
line defines the update matrix Ty = — (i7y — S^y). 


The rules for evolving the psips, which resemble those 
used in FCIQMC^^, follow from Eq. (5): 

1. Psips can spawn from a density matrix element pik 
to py with probability Ps(ik —>■ ij) = A/3|rkj|/2, 
with sign(py) = sign(pik) x sign(Tkj); a similar 
spawning process takes place from pkj to py. 

2. Psips on the density matrix element py clone/die, 
whereby their population is increased or decreased, 
with probability Pd(ij) = ^|7ii + ^jjl- The popu¬ 
lation is increased if sign(Tii + Ijj) x sign(py) > 0 
and decreased otherwise. 

Additionally, we annihilate psips of opposite sign on 
the same density matrix element to improve the effi¬ 
ciency of the algorithm and overcome the Fermion sign 
problem^^’^^. We note that, unlike PIMC methods where 
the quality of averages depends on the average sign of 
the sampled paths®’in FCIQMC and DMQMC, we re¬ 
quire a system specific and basis set dependent critical 
psip population to obtain correct low temperature and 
ground state estimates^^’^^^^®. 

The simplest starting point for a simulation is at /3 = 0, 
where the density matrix is the identity and can be sam¬ 
pled by occupying diagonal density matrix elements with 
uniform probability. A simulation then consists of prop¬ 
agating the initial distribution of psips with the rules de¬ 
scribed above to a desired value of /3. Estimates for ther¬ 
modynamic quantities can be found by averaging over 
many such simulations, a single one of which we call a 
‘/3-loop’. 

In Fig. 1(a) we see that a direct application of 
DMQMC to the dense UEG can result in estimates for 
the internal or total energy that are too high in the inter¬ 
mediate temperature range. This can be understood by 
noting that, at = 1, the ground state of a few-electron 
UEG system is well described by a single (Hartree-Fock) 
determinant, |Do). The probability of initially selecting 

this determinant, however, is (^) , which rapidly ap¬ 

proaches zero as M increases. If, by chance, the Hartree- 
Fock determinant or another low-energy determinant is 
sampled at /3 = 0, the population of psips arising from 
that low-energy determinant will dominate the simula¬ 
tion, but most simulations miss the low-energy part of the 
Hilbert space altogether. As shown in Fig. 1(b), this sam¬ 
pling problem reduces as the number of /3-loops (or the 
population of psips per /3-loop) increases, thus increasing 
the chance of sampling the low-energy space; however, 
this approach soon becomes impractical. Fig. 1(a) shows 
that, by moving to the interaction picture, we can effec¬ 
tively solve this sampling issue and regain FCTquality 
thermodynamic averages. 

B. Moving to the Interaction Picture 

There are two sampling issues present when treating 
real systems; the distribution of weight in the density 
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FIG. 1. Panel (a) shows the total energy calculated for a 
seven-electron spin-polarized electron gas at ra = 1 with 
M — 33 plane waves in the total momentum K = 0 subspace 
using DMQMC and an initial psip population Np = 10"^. In¬ 
creasing the number of /?-loops, Np, from 10^ (squares) to 10^ 
(crosses) results in a more accurate answer being reproduced. 
We also see that the error bars do not reflect the true errors for 
Np = 10^ in the intermediate /3 regime. Panel (b) shows that 
the average occupation on the Hartree-Fock density matrix 
element (|IIo)(I3o|) is under-represented in the intermediate 
temperature range. Also plotted in (a) are the IP-DMQMC 
results (circles) for the total energy calculated using Np = 10^ 
and only 10 /3-loops (see Section II B). 


where we have used the usual definition of an operator 
in the interaction picture: 

Ai{t) = (11) 

In practice the exponential factors appearing in Eq. (10) 
are time-consuming to evaluate and we prefer to work 
with Eq. (9). Since we choose if® to be diagonal in our 
determinantal basis set of Hartree-Eock eigenstates, the 
only modification to the original DMQMC algorithm is 
that Pd(ij) = — Hjj\. This scheme, which we 

dub interaction-picture DMQMC (IP-DMQMC), has the 
added benefit that there is typically no death along the 
diagonal as long as > Ha. This overcomes a related 
issue in DMQMC simulations of large systems, whereby 
the weight of walkers on the diagonal decays nearly to 
zero with /3; this was previously remedied with the use 
of importance sampling^^. 

Eq. (8) shows that IP-DMQMC only samples the cor¬ 
rect distribution at r = /3 so that, unlike in DMQMC, 
where the whole temperature range is sampled in one 
simulation, separate simulations are required for each f3 
value. As in DMQMC, estimates for observables require 
averaging over multiple independent simulations. 

Referring back to Fig. 1 we see that working in the in¬ 
teraction picture effectively eliminates this sampling is¬ 
sue, with the correct total energy being reproduced using 
small numbers of psips and /3-loops. 


matrix changes rapidly as a function of /3, and important 
determinants are rarely present in our initial configura¬ 
tions. Feynman originally pointed out that if we can 
write H = + V, where V is small compared to , 

then the quantity p will be a slowly-varying function 
of /3^^. This does not solve the issue of selecting impor¬ 
tant determinants, so we define an auxiliary matrix 


/(r) = 

(6) 

which has the properties 


f^r = 0) = e-^^\ 

(7) 

f{r = /3) = = m- 

(8) 


From Eq. (7) above we see that, by working with the op¬ 
erator /, we can start the simulation from instead 

of the identity. For most weakly-correlated systems this 
should provide a good first approximation to the distri¬ 
bution of weight in the fully interacting density matrix. 

Differentiating Eq. (6) with respect to r we find 


= H°f-fH 

(9) 

II 

1 

1 

(10) 


C. Sampling the initial condition 

The choice of is somewhat arbitrary, but it should 
allow for an efficient sampling of /(t = 0) and this is 
most easily achieved if is non-interacting. In princi¬ 
ple, any initial density matrix can be sampled using the 
Metropolis algorithm^®, but we have found this approach 
problematic due to the long equilibration times required 
at low temperatures and in large basis sets. An alter¬ 
native method, which is free from such issues, is to use 
knowledge of the grand canonical density matrix corre¬ 
sponding to and sample this in such a way that the 
desired, canonical, distribution is reached. 

Consider the grand canonical density matrix 


Pgc — c 


-P{Ho-pN) 


( 12 ) 


where Hq = is some non-interacting Hamilto¬ 

nian whose single-particle eigenvalues Si are known, N 
is the number operator and fi is the chemical potential, 
which can be determined from the implicit relationship 

{N)gc = E + 1 


(14) 
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where we have identified the usual Fermi-factor, pi. The 
grand canonical partition function, Zgc^ can be written 
as 


N {nf} 


where {n^} denotes a set of occupation numbers such 
that ^’^cl Hi € {0,1} for fermions. 

The probability of selecting a particular set {nf} is 


PGc{{n^}) 


1 

Zgc 




(16) 
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However, we wish to generate determinants in the canon¬ 
ical ensemble where the correct probability is 

Pc{{nn) = ^ n (17) 

and 

Zc = E (1^) 

{"fi 

We see that Pciin^}) oc PGcilnf^})- Thus, by indepen¬ 
dently occupying orbitals with probability pi and then 
discarding those configurations with {N) ^ N, we attain 
the correct proportionality factor ZggI'Z'C- Only about 
one in '/N of the configurations sampled has the right 
value of IV, but the sampling process is so fast that very 
little computer time is required for any system size in 
the reach of many-body simulation methods. The chem¬ 
ical potential can be obtained by numerically inverting 
Eq. (13) in the appropriate finite basis set, a procedure 
we carry out using Scipy^®. A demonstration of the above 
procedure is given in Fig. 2, where we see that (iJ) is in¬ 
deed a slowly varying function of t and that the correct 
estimate is reproduced at r = /?. 

Finally, we note that any diagonal density matrix can 
be obtained by reweighting the configurations which re¬ 
sult from the above sampling procedure as 

PneAin^}) = J^oid({nf (19) 

where E^ew and Eoid are the new and old total energies 
of a given configuration {nf^}, respectively. 


III. BASIS SET EXTRAPOLATION 

To treat the UEG using DMQMC we need to work in a 
finite basis set of M plane waves; thermodynamic quanti¬ 
ties will therefore need to be extrapolated to the M ^ oo 
limit. At T = 0 it has been found^®’^^’^® that the cor¬ 
relation energy, = E — E^p, for unpolarized systems 
converges like M~^. Recent CPIMC results^® were also 


FIG. 2. Variation of {H) with r using 77° = X^k^kC^Cki i-®-: 
the free-electron expression with ep — |k^. The grand canon¬ 
ical procedure described in Section IIC was used. The system 
shown is a four-electron, spin polarized gas at = 1, M = 81 
and d = 1- For these results we used approximately 10° psips 
and averaged over 100 simulations. The dashed line repre¬ 
sents the exact FCI result, which IP-DMQMC reproduces at 
T = (3 as expected. For comparison (77) = 50.751(4) Ha at 
13 = 0. 

obtained using an M~^ extrapolation, although in prin¬ 
ciple this relationship only holds for an unpolarized sys¬ 
tem; we have found that, for polarized systems, extrap¬ 
olating with M”®/® results in a better fit®°, consistent 
with similar observations by other authors®°’®^. Based 
on the analysis presented here, Schoof et al.^^ have used 
the extrapolation with CPIMC. In addition, we 

find that the convergence of the total energy is strongly 
dependent on temperature. 

At T > 0 there is a competition between the conver¬ 
gence of the kinetic and potential energies with M. To 
investigate this further we focus on a two-electron spin- 
polarized system, which can be solved exactly using diag- 
onalization in large basis sets. In Fig. 3 we see this com¬ 
petition between energy scales: the total energy initially 
increases rapidly with basis-set size before appearing to 
saturate. As the size of the basis set is further increased, 
a slight reduction in the total energy is observed, with 
the residual error apparently proportional to M”®/® (see 
inset of Fig. 3). In this regime the convergence of the 
total energy is dominated by exchange and correlation 
effects. 

The initial increase of the total energy with respect 
to M at non-zero temperatures can be understood by 
looking at the non-interacting total energy as a function 
of basis-set size, which is most easily analyzed in the 
grand canonical ensemble. The non-interacting basis-set 
error is 

AEo{M) =Eoioo)- EoiM) (20) 

= E (21) 

k>kc 
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FIG. 3. Behavior of the FCI total energy with basis-set size 
for an N = 2, Vs = 1 spin-polarized system at 0 = 0.5. Here 
we see the competition between the exponential convergence 
of the total energy at low M (main plot) and the 
behavior for high M (inset). 


where kc is the plane-wave cutoff. For » 1, 

this can be approximated as 

pOO 

AEo^ / (22) 

Jec 

where we have used pk « for ^ M- Hence 

poo 

AEq{M)k / + (23) 

Jo 

poo 

«e3/2e-/3(ec-M) / e-^^dx, (24) 

Jo 

if £c ^ so that (Ec + x)'^/'^ « everywhere 
is significant. It then follows that the leading-order cor¬ 
rection is 


AF;o(M) « (25) 

From Eq. (25) we see that the kinetic energy begins 
to converge exponentially once Ec ~ jd~^ or M « 
« iV03/2, where V is the simulation-cell vol¬ 
ume. In practice, we find that for large 0 the kinetic 
energy and hence the total energy converge quite slowly. 
This is an issue for DMQMC simulations as the cost of a 
calculation increases with basis-set size. 

We can mitigate some of these issues by instead ex¬ 
trapolating the temperature-dependent correlation en¬ 
ergy, Ec{P,M) = E{P,M) - EiiF{l3,M). The infi¬ 
nite basis-set total energy can then be reconstructed as 
E{/3, M = oo) = Ehf(/3, M = oo) -|- Ec{(3, M = oo). We 
calculate the ‘Hartree-Fock’ energy, Ehf = (77)hf, using 
the density matrix^^ 


PHF = 5]e-''^"|A)(A|, (26) 


FIG. 4. Exponential convergence of the Hartree-Fock total 
energy with basis-set size for A = 4, = 1 and 0 = 4. Note 

that no Madelung contribution is included for the Hartree- 
Fock estimates in this figure. The infinite basis-set limit for 
i?HF is estimated as 70.792(1) Ha. For comparison, we plot 
{H)o and {Ho)o, where the trace is now with respect to the 
non-interacting density matrix. Also plotted is the total en¬ 
ergy calculated in the grand canonical ensemble using a finite 
basis set, i.e., {Hq)gc = lOk'" £kPk. 


where Ef^^ = {Di\H\Di) and the sum runs over all de¬ 
terminants in the basis set. £’hf(/ 3) can be found using 
the sampling procedure outlined in Section IIC, i.e., 

^HF . 

1 

= y^£:.HFg-/3(iSr"’-B?)g/3B[’ 

^HF . 

1 

SiW(i)p(i) 

where u;(i) = and p(i) = . Thus, 

by generating determinants as described in Section IIC 
and reweighting them using r/;(i), we can instead sample 
Phf and, as a result, estimate Ehf as desired. In Fig. 4 
we show the convergence of E-hf{P,M) as a function of 
basis set for a four-electron, spin-polarized system at = 
1 and 0 = 4. Note the large basis-set sizes required 
to converge the total energy to within statistical error 
bars. Fig. 4 also shows various other ‘non-interacting’ 
or ‘mean-field’ energy estimates as functions of M. Any 
of these could in principle be subtracted from E{f},M) 
to define a correlation energy, but the quantity defined 
by subtracting i?HF(/?, M) extrapolates most smoothly to 
the infinite M limit. The non-interacting grand canonical 
energy, {Ho)gc, is significantly larger than the canonical 
estimates. 

Fig. 5 shows how A(/3) M) depends on M at a number 
of different temperatures. For small basis sets A shows a 
power-law decay with M, but this ceases for large enough 
M and the energies begin to increase again. We believe 
that the increase is due to kinetic effects that are not 


(27) 

(28) 
(29) 
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FIG. 5. Comparison of the convergence of Ec with basis-set 
size calculated using exact diagonalization for a two-electron 
spin-polarized system at rs = 1 for different values of 0. The 
0 = 2 data has been shifted down by 0.0002 Ha for visibility. 
The dependence of Ec on 0 is non-monotonic. The converged- 
basis-set correlation energies at 0 = 8 and 0 = 2 are very 
similar, but the correlation energy at 0 = 4 is more negative. 

captured in the non-interacting expression we subtract. 
The value of M at which the correlation energy begins 
to increase again corresponds to the onset of the power- 
law convergence of the total energy observed in Fig. 3. 
This non-variational behavior of the internal energy with 
respect to M is not surprising as, at finite temperatures, 
it is the free energy that satisfies a variational principle. 

We can estimate Ec{P,M = oo) by taking the value 
calculated with the largest basis set after the minimum is 
reached. This will in general over estimate Ec (it will be 
too negative) but the remaining discrepancy is typically 
smaller than the stochastic error bar. The systematic 
errors left after extrapolating Ec{l3,M) in this manner 
are well within chemical accuracy (~ ksT) and can be 
orders of magnitude smaller than those introduced by a 
direct extrapolation of the total energy. 


IV. RESULTS 
Computational Methods 

All calculations discussed in this paper were performed 
using the HANDE code^^. Unless otherwise stated, the 
QMC calculations used real amplitudes^^^^’^ to sample 
the density matrix, which improves stochastic efficiency 
compared to integer weights. The full data set is available 
in the supplementary material^®. 


A. Four electron uniform electron gas 

Using the procedures outlined above we are now in 
a position to provide exact benchmarks for the UEG 



P (Ha-i) 

FIG. 6. Comparison of IP-DMQMC results (markers) with 
FCI results (dashed-lines) for the UEG with A = 4 and Cs = 1 
in two different basis-sets. The inset shows the low T behavior 
where we see increasing the basis-set size serves to decrease 
the total energy in contrast to the high T behavior, where the 
opposite occurs. These calculations used integer rather than 
real weights and were run using approximately 1000 psips and 
we averaged over 100 simulations. 

in small simulation cells across the relevant parameter 
space. In this first study we focus on the four-electron 
spin-polarized system, which is the smallest non-trivial 
system and one for which there already exist bench¬ 
mark calculations^®. All energies contain the Madelung 
contribution®® where appropriate. As a first step we com¬ 
pare our four-electron IP-DMQMC results to FCI results 
in small basis sets and see perfect agreement across the 
whole temperature range (Fig. 6). 

We have extended these results to basis sets far beyond 
the reach of conventional full diagonalization procedures; 
the largest space sampled here contains approximately 
10®® density matrix elements. We used the initialization 
procedure outlined in Section IIC and the free-electron 
Hamiltonian for H® for rg < 1; for rg > 1 we found it ad¬ 
vantageous to use Hartree-Fock density matrix defined in 
Eq. (26). The calculations were initialized with 10®-10^ 
psips and the results averaged over 100-5000 simulations, 
each using a different random number seed. Time steps 
At ranging from 0.01/Ep to 0.001/Ep were used, with 
a smaller time-step required at lower Tg; the values cho¬ 
sen were small enough that we could resolve no time-step 
error within the statistical errors. Each (rg,0,M) cal¬ 
culation was typically run for 2 hours on 48 cores with 
a total computational cost of approximately 80000 core 
hours. The separate calculations of Uhf required 9000 
core hours. 

In Figs. 7 and 8 we show the convergence of the IP- 
DMQMC results with basis set at low and high temper¬ 
atures, respectively. We note that the behavior found in 
the two-electron system is also found in the four-electron 
system; in particular the non-trivial dependence of the 
correlation energy on M is reproduced. We find that a 
direct extrapolation of the total energy with respect to 
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M is best for 0 < 0.25 as it is here where kinetic ef¬ 
fects are minimal and there is a clear trend in the total 
energy for the basis set sizes considered. The procedure 
outlined in Section III is best suited for temperatures 
above this, becoming increasingly useful for 0 > 2 as 
more highly excited states become accessible, requiring 
prohibitively large basis sets for a direct extrapolation to 
be possible. In between these too regimes both methods 
produce statistically identical results^®. Fig. 9 summa¬ 
rizes our results and shows perfect agreement with the 
available CPIMC data from Ref. 29. Further results at 
higher temperatures and other rg values are available in 
tabular form in the supplementary material and again 
agree with the available CPIMC results. 
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FIG. 7. Total energy of a four-electron spin-polarized sys¬ 
tem at rs = 1 and 0 = 0.0625, showing a convergence with 
M-5/3. The dashed line represents an extrapolation to the 
infinite-basis-set limit carried out using a least-squares fit 
as implemented in Scipy^®. At this low temperature, a di¬ 
rect extrapolation of the total energy works better than the 
correlation-energy extrapolation technique discussed in Sec¬ 
tion III. 


FIG. 8. Convergence of Ec with basis-set size for the four- 
electron system at Cg = 5 calculated using IP-DMQMC. Note 
the similar behaviour to that found in the two-electron case 
in Fig. 5. At the high temperatures considered here, the 
correlation-energy extrapolation technique introduced in Sec¬ 
tion III works much better than the total-energy extrapola¬ 
tion illustrated in Fig. 7. 
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V. DISCUSSION AND CONCLUSIONS 

In this paper we have demonstrated how DMQMC can 
be applied to realistic systems. By moving to the inter¬ 
action picture we have removed sampling issues found 
when treating weakly-correlated systems with large basis 
sets. 

We have examined in detail the convergence of the to¬ 
tal and correlation energies with respect to basis-set size 
M and temperature using a system accessible to FCI cal¬ 
culations. We found that, in general, these quantities ex¬ 
hibit a non-trivial dependence on M attributable to the 
competing energy scales present. By developing a simple 
Monte Carlo sampling scheme, we showed that it is possi¬ 
ble to reduce the error in extrapolating these quantities to 
the complete basis-set limit by at least an order of magni¬ 
tude at high temperatures. We believe this analysis and 


FIG. 9. (a) Extrapolated total energies per particle for the 
four-electron system at = 1,2,4 showing exact agree¬ 
ment with the GPIMG results of Ref. 29. Dashed lines 
are meant as guides to the eye. (b) Relative deviation, 
(Edmqmc — Ecpimc)/Edmqmc, as a function of tempera¬ 
ture showing statistically identical results for Vs = 1 (circles). 
Vs = 2 (diamonds) and rg = 4 (squares). Further results at 
higher temperatures and other Cg values are available in tab¬ 
ular form in the supplementary material and again agree with 
the available GPIMG results. 


developments will be useful when treating molecular sys¬ 
tems at finite temperatures. In addition, our approach 
to calculating the ‘Hartree-Fock’ energy should a useful 
first approximation when providing accurate benchmark 
calculations for systems away from the thermodynamic 
limit and in analyzing single-particle finite-size effects for 






































the UEG at non-zero temperatures. 

Using these developments we have reproduced the four- 
electron CPIMC benchmarks of Ref. 29 and provided re¬ 
sults at higher temperatures. We hope that these small- 
system results will aid the analysis of the apparent dis¬ 
crepancies between other QMC methods for larger sys¬ 
tem sizes^° and serve as benchmarks for other QMC 
methods based in configuration space. 

Whilst the results presented here are for much smaller 
systems than those accessible by RPMIC and CPIMC, 
DMQMC provides access to exact finite-temperature 
data for a given basis set. The main limitation on the 
system size is the critical population (determined by the 
plateau height^^’^^) required to sample the density ma¬ 
trix. There are several grounds for optimism. The sign 
problem is much weaker at higher temperatures, imply¬ 
ing that larger systems will be accessible albeit over a 
restricted temperature range. Further, in our previous 
study^^ we found that the plateau height in DMQMC, 
which is a measure of the strength of the sign problem, 
was roughly the square of that in FCIQMC, which would 
suggest a rather limited utility of our method. Remark¬ 
ably, however, we have found that, for the UEC at var¬ 
ious Ts values and for some other models, the plateau 
height is only a small multiple of the FCIQMC plateau 
height. For example, for = 4, Ts = 5 and M = 1045, 
the FCIQMC plateau height is roughly 8000 psips in com¬ 
parison to 90000 in IP-DMQMC at 0 = 0.0625. Finally, 
there is evidence that methodological developments will 
increase the scope of systems that can be treated with 
DMQMC. Civen that the initiator approximation^^ en¬ 
abled FCIQMC to be applied successfully to large UFC 
systems across a range of densities^®’^^, we are confident 
that, using similar approximations and developments in 
importance sampling, DMQMC will become a competi¬ 
tive method in treating degenerate Fermi systems. 
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